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Advisor: Dr. Farzad Mashayek 



A detailed numerical study of evaporating capillary jets is presented. The analysis is 
performed through use of a Galerkin finite element method with penalty formulation for solving 
the equations of motion and a flux method for tracking the free surface. 

A parametric study is performed to analyze the temporal instability of the evaporating jet. 
Through varying the evaporation rate, Reynolds number, disturbance wavenumber, initial 
disturbance amplitude, and density ratio the outcomes of jet breakup are investigated. Also, 
pressure distribution inside the jet and multiple satellite drop formations are analyzed. Results 
are compared to existing analytical conclusions made from linear stability analysis. 

This study reveals that surface evaporation has a destabilizing effect for the low-speed 
jets, which are considered here. That is, evaporation flux is greater at the neck than the crest, 
which accelerates the wave growth. Satellite drops also reduce in size as evaporation rate is 
increased. This reduction is seen in both the radial direction due to vapor leaving the surface and 
along the axis of symmetry due to decreased breakup time. 
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Nomenclature 



Q 

h 

I 

k 

K 

I 
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m 

n 

N 

P 

Pr 

Q 

r 

R 

Re 
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t 

T 

T 

M 

u 

V 

V 

V 

w 

We 

z 



specific heat 
height function 
identity matrix 
wavenumber 
jet curvature 
length 

latent heat of vaporization 
evaporation flux 
outward unit normal 
shape factor 
pressure 
Prandtl number 
evaporation rate/length 
radial coordinate 
instantaneous radius 
Reynolds number 
density ratio 
time 

temperature 

Newtonian stress tensor 
jet axial velocity 
free stream velocity 
radial velocity 
velocity vector 
volume 

velocity vector 
Weber number 
axial coordinate 



z 


Ohnesorge number 


p 


dimensionless evap. rate 


r 


surface area of element 


5 


change in 


e 


disturbance amplitude 


T| 


surface wave height 


K 


thermal conductivity 


X 


wavelength 


V 


kinematic viscosity 


P 


density 


o 


surface tension coefficient 


X 


stress tensor 


0) 


growth rate 


Q. 


volume of element 


Subscripts 

b 


boiling 


c 


curvature 


e 


element 


8 


gas 


ij 


node number 


1 


liquid 


n 


neck 


s 


swell, surface 


z 


derivative 


0 


initial value 



V 



1. Introduction 



The instability of liquid jets has been a classical fluid mechanics problem for more than a 
century. Early investigations focused primarily on understanding the mechanisms of jet breakup 
while recent jet instability studies have concentrated on more accurate predictions of growth rate, 
drop size after breakup, and effects due to variations in initial disturbance amplitude, disturbance 
type, and fluid properties. Despite extensive previous studies cited in the literature and briefly 
reviewed below, it has been only recently that the dynamics of a jet undergoing surface 
evaporation has been investigated by Lian and Reitz (1993) using a linear stability analysis. 
While this study reveals several new and interesting phenomena, it is, understandably, unable to 
predict accurate results near the breakup point due to the assumption of infinitesimal perturbation 
of the liquid jet. 

The objective of the present work is to relax some of the assumptions adopted by Lian 
and Reitz (1993) and to provide a more versatile computational analysis of evaporating jets by 
implementing a Galerkin finite element method in conjunction with a height flux method (HFM) 
for interface tracking, which was developed by Mashayek and Ashgriz (1993). The 
dimensionless continuity and momentum equations are solved on a moving mesh consisting of 
four-node quadrilateral elements. Considering an axisymmetric incompressible Newtonian 
liquid jet of infinite length initially subjected to a periodic surface disturbance a parametric study 
is performed to show how evaporation affects capillary jet instability while varying liquid 
Reynolds number {Re), disturbance wavenumber {k), disturbance amplitude (£o), evaporation rate 
(P), and density ratio (s). Similar to Lian and Reitz we focus on the temporal evolution of the jet 
thus only one wavelength is examined at a time - in practice, we solve for only one half of this 
wavelength due to the symmetry of the jet. Also similar to Lian and Reitz (1993), the surface 
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evaporation of the perturbed jet is approximated by that of a spherical drop with the identical 
curvature as the local surface curvature of the jet. In this manner, there is no need to solve for 



the gas phase. 

The observation that a liquid jet issuing from a nozzle will eventually breakup into small 
drops when subjected to even minute disturbances led Rayleigh (1879) to provide the first 
analytical description for the temporal instability of an inviscid, incompressible jet using a 
normal mode analysis. He showed that an axisymmetric harmonic disturbance of the form: 



where (0 is the growth rate, t is the time, and 7o and I\ are the modified Bessel functions of the 
first kind. In (1) and (2), length and time variables have been non-dimensionalized by Vo and 
To/ve, where Ve is the capillary-wave velocity (o/pro)'^^ Here, ro is the undisturbed radius of the 
jet, a is the liquid surface tension coefficient, and p is the density of the liquid. This result 
predicts that disturbances grow in time for k<l and oscillate for /:>!, and the maximum growth 
rate occurs at k=0.691. In other words, if the wavelength (h=2KrJk) of an axisymmetrical 
disturbance is greater than the cylinder circumference the jet will be unstable and if the 
disturbance wavelength is less than the cylinder circumference the jet will be stable. 

For over a century since Rayleigh noticed that surface tension had to work against inertia 
in causing instability researches have continued to extend his work in both theoretical and 
experimental realms. For example, viscous effects have been considered by Weber (1931), who 
analyzed the motion of thin sections of fluids, and Chandrasekhar (1961), whose book focused 



ri = e„ exp((ot - ikz). 



( 1 ) 



would grow in time according to: 




( 2 ) 
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solely on the subject of instability using the Navier-Stokes equation to examine it. To test the 
results of linear theory quantitative experiments were performed by such researchers as Goedde 
and Yuen (1970). The experiments clearly identified that drop size varied at breakup leading to 
the conclusion that nonlinear behavior became dominant during instability. Yuen (1968) was the 
first to extend Rayleigh’s linear analysis by considering a third-order perturbation expansion, 
which would give some explanation of the production of satellite drops between the larger main 
drops. Further experimental investigation of both satellite and main drops have been done by 
Vasallo and Ashgriz (1991) and Kowalewski (1996), which looked at the production of multiple 
satellites and the description of the jet surface near breakup, respectively. 

Due to recent development of computational facilities and numerical techniques flow 
simulations have become common place. Inviscid flow was studied by Mansour and Lundgren 
(1990) and Stokes flow has been considered by Stone and Leal (1989). These two simplifying 
factors breakdown close to breakup because both viscous and inertial effects become important. 
In this case the full Navier-Stokes equation has to be considered in the fluid domain. A finite 
difference one-dimensional model was considered by Eggers and Dupont (1994) to simulate drop 
formation. The most extensive study of jet breakup using the Navier-Stokes equation was 
performed by Ashgriz and Mashayek (1995) using a Galerkin finite element method. For a 
review of the computational methods used for free surface flows see Tsai and Yue (1996). 
Further reviews of the problem of jet instability are given by Bogy (1979) and more recently by 
Eggers (1997). 

Unfortunately, little mention is made about evaporating jets when studying the problem 
as a whole. This is primarily due to the complexity of determining an evaporation equation for 



3 









i-f . 












*»? ^ J 



A 



% 



# '* '*■ ■•■. 









a>' — *4'^ 



♦ 



.J 



• I? s 1 » ' Vic*' 

11*-^ < ^ - • I . ' MV I ’ I , I’ 



IMP i 



* ^■■. -r-r v:^ 



tt. » 



f <1 «■ 1 M 



mnHitfui |ii. 4l«*m aMf«o mu ffwn^w^w^ r 







H'. '''*' js^' 



ir 



K5* ..^ .. 



an irregular shaped surface. Hopefully future experimentation will be performed to verify the 
spherical model used herein. 



2. Formulation and Methodology 

Consider a viscous liquid circular jet injected with velocity U into an inviscid gas at 
temperature in zero gravity. Since we are dealing with interactions of the gas with only one 
jet, it can be assumed that the total heat capacity of the gas is much larger than that of the liquid 
and TL,, far from the surface of the jet, is constant. The implementation of an evaporation model, 
described below, eliminates the need for the solution of the energy equation; thus the governing 
equations include only mass and momentum conservations. Assuming both the liquid and the 
gas are incompressible with constant properties, these equations are described as: 

Vw, =0, (3) 

^^ + w. -Vw,. = — ^-Vp, +V;W,., i = l,g (4) 

dt p,. 

where / and g stand for liquid and gas, respectively. In (3) and (4), w, p, p, and v denote the 
velocity, pressure, density, and kinematic viscosity, respectively. For a liquid jet evaporating 
with a flux m , the jump condition and the normal stress balance on the interface are given by 
(Lian and Reitz, 1993): 



m = p,(w, -wj n = p^(w^ -wj n. 


(5) 


m(w^ - w, )• n + [n • (t^ - 1, )J n + aV • n = 0, 


(6) 



where a, t, Wj, and n denote the surface tension coefficient, the stress tensor including the 
pressure component, the interface velocity, and the outward unit normal, respectively. The first 
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term on the left-hand side of (6) is the recoil force on the jet due to surface evaporation. This 



force stems from a net momentum that is impacting the surface. Eliminating Wj in (5) yields: 



pressure may be neglected in comparison to the effects of the pressure fluctuations in the liquid. 
Further, the mean pressure of the gas affects only the mean pressure of the jet and, for 
incompressible liquid, does not influence the instability of the jet. As a result, we can neglect the 
ambient gas pressure and, for inviscid gas, eliminate Xg from (6). Then, by implementing (7), (6) 
yields: 



This simplifies the problem significantly, as there will be no need to solve for the gas phase. It 
must be noted that, in this manner the mean pressure of the jet is calculated relative to the mean 
pressure of the ambient gas. 

There is no analytical solution for evaporation of a jet with a perturbed surface. 
However, following Lian and Reitz (1993) one may approximate the unit area evaporation rate 
of a wavy curved surface of local radius of curvature Rc with that of a sphere of the same radius. 
In this study, assuming that the surface temperature of the jet is always at the boiling temperature 
Tb, the flux of evaporation is modeled as: 




(7) 



In this study we consider Pg«p/, therefore the effects of the fluctuations in the ambient gas 



_ m 

n • t, • n = aV • n H . 



( 8 ) 



OT = -^ for R^>0, 



(9) 



where. 
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(l + 0.276 



where K, Cp, and denote the thermal conductivity, the specific heat, and the latent heat of 
evaporation, respectively. The Reynolds number. Re, and the Prandtl number, Pr, are defined as 
Reg-lUrJVg and 7*r^=p^VgCp^/Kg. Equations (9) and (10) are a modified form of the J^-law for 
evaporation of a single drop in zero gravity (Williams 1985 and Turns 1996). The modification, 
here, is by substituting the local radius of curvature of the perturbed jet for the radius of the 
spherical drop. This modification was first introduced by Lian and Reitz (1993). The physical 
reasoning is that the surface of the deformed jet may be “locally” considered as the surface of a 
spherical drop having the same radius of curvature. It is noted that the utilization of (9) is 
limited to surface deformations such that Rc > 0. 

In the remaining equations (for the liquid phase) an axisymmetric incompressible 
Newtonian liquid jet in zero gravity is considered. The variables are non-dimensionalized by the 
undisturbed jet radius r^, and the characteristic time {pro /of The dimensionless continuity and 



momentum equations are: 



Re 



Vv = 0, 

3v 

^ + vVv |=V-T, 
3t 



( 11 ) 



( 12 ) 



where v=W/=(m,v) is the velocity vector and T=t/=-pI+[Vv+(Vv)^] is the stress tensor for a 
Newtonian fluid with p now denoting the normalized liquid pressure. The Reynolds number in 
(12) is defined as Re={\Ni){<5rJpi)^^. Non-dimensionalizing (8) the stress balance on the 
interface yields: 

Tn = Re{K + aK^)n, (13) 
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where K is the curvature of the surface and 



oc = —, p = ^ 



s = ■ 



Pi 



( 14 ) 



As a result of the non-dimesionalization, the parameters affecting the instability of an 
evaporating jet reduce to Reynolds number Re, the dimensionless evaporation rate P, and the 
density ratio s. In the next section we investigate the effects of these parameters as well as the 
effects of initial disturbance amplitude £o and wavenumber k via numerical simulation. 

In this investigation it is assumed that the free surface can be represented by a height 
function h{z,t), as shown in Fig. I, thus the curvature K is calculated using: 



K = 



1 



( 15 ) 



Since a temporal analysis is considered symmetry boundary conditions can be applied on the axis 
and planes of symmetry: 




Figure 1 : Ruid subvolumes represented by their height and thickness. 
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— = 0 and V = 0, at r = 0 
or 



( 16 ) 



u=0 and — = 0, at z = 0,— . (17) 

dz 2 

A Galerkin finite element method is used to solve (11) and (12). Using a penalty function 
formulation (Hughes, Liu and Brooks 1979), the pressure term is eliminated from the set of 
unknown variables by absorbing the continuity equation into the momentum equation. In this 
formulation the pressure is defined as: 

p = -YVv, (18) 

where Y is a large number of order (10^) depending on viscosity and Reynolds number. Four- 
node bilinear isoparametric elements are used to approximate tbe velocity distribution over each 
element: 



v(r,z,U=5^V((UA,.(r,z,r). 



(19) 



1=1 



A moving mesh is considered to discretize the computational domain making the shape 
functions, time dependent. In order to obtain the finite element formulation, (12) is 
multiplied by tbe shape function Nj and integration is carried over the element volume. After the 
divergence theorem is invoked the resulting equation is: 



f NjRe 
Jowl ■' 






^+v Vv +VA^J •[-pI + ^v + (Vv)^| |(fQ = NjT-ndT, 
dt j ^ fw 



( 20 ) 



where Q is the volume and T is the surface area of the element. Through substitution of (13) and 
(18) into (20) we obtain the following closed form finite element formulation: 



f \N^Re 
Jow ■' 



3v 

— + V- Vv 
dt 



\ 



+ VN] • (V • v)l + [s7v + (Vv/ ] ](/Q = N.Re{K + aK^)\dr. (21) 
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The above formulation is based on a fixed mesh. In the present problem the mesh is moving, 
therefore special treatment of the time derivatives is required. Since the shape function is time 
dependent, the time derivative of velocity in discretized form becomes: 




( 22 ) 



The last term of (22) introduces a new convective term in (21). Here, we allow the motion of the 
nodes in the r-direction only, according to the following simple rule: 



where the subscript i refers to the node number, and c=c(z,0 is a constant for each column of 
nodes in the radial direction defined as: 



In this case, Mashayek and Ashgriz (1993) have shown that the total derivative of velocity 
becomes: 



Substitution of (25) into (21) completes the finite-element formulation for the present problem. 

The free surface of the jet is determined using the Height Flux Method (HFM) developed 



subvolumes of width 5z/ and volume V/. The location of the free surface on the left and right 
hand sides of a subvolume is given by hi and /i,+i at time t, respectively. At the end of each time 
step fluxes of the fluid from each subvolume to its neighboring subvolumes are calculated using 
the velocity field determined by the finite element solution of the governing equations. Next, 
with the knowledge of the evaporation rate, the volume change 6K of each subvolume due to 
evaporation at the interface is calculated using: 



z,(r-H50=Z/(0=constant, ri(t+bt)=cri(t). 



(23) 



c = h{z,t + 8t)lh{z,t). 



(24) 




(25) 



by Mashayek and Ashgriz (1993). As shown in Fig. 1 the domain is divided into several vertical 
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h.^, +/j, \ 



^ |/irp6/6/, 
^ ) 



(26) 



where 5/ indicates the length of the interface confined within the subvolume. The volume of 
fluid within each subvolume is further modified by calculating the flow through the side planes. 
Then, it is assumed that the part of the interface located between any neighboring pair of 
subvolumes can be approximated by a line segment described as h=az+b, where a and b are 
constants which can be determined by the known volumes. More details and the accuracy of this 
method are discussed by Mashayek and Ashgriz (1993). 



Mesh Size 


4x30 


4x40 


4x50 


4x60 


Growth Rate 


0.712 


0.712 


0.712 


0.712 



Table T. Typical convergence test for growth rate: Re=200, £o=0.05, k=0J, P=3xl0‘^, and 5=10’^. 



Before numerical results were compiled a convergence study was performed to determine 
the effect of mesh size on growth rate. The results are shown in table 1. From this we conclude 
that grid size does not affect growth rate since only the height of the free surface at the neck and 
swell where of importance for calculating this value. To save computational time the majority of 
cases were ran with a minimal number of elements. The mesh size for the determination of 
growth rate was limited to four elements in the radial direction and the axial elements were 
changed based on the wavenumber as given in table 2. For the cases where pressure distribution 
and velocity vectors were plotted the grid was refined in order to obtain better print quality and 
color dispersion. 



Wavenumber 


0.3 


0.4 


0.5 


0.6 


0.7 


0.8 


0.9 


Number of elements 


50 


40 


40 


30 


30 


20 


20 



Table 2: Number of axial elements used in simulations for determination of growth rate. 
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3. Effect of Evaporation on Jet Instability 



We consider the temporal evolution of a liquid jet, initially at rest with a spatially 
harmonic surface described as: 

r=R-eoCos(kz), (27) 

where R is determined such that the volume of the jet is kept constant when the initial amplitude 
is changed: 



R = (l-ie;f. (28) 

Due to symmetry, only half a wavelength is considered. The trough (neck) of the initial surface 
is set at z=0 and the crest (swell) at z=X/2. The dynamics of this jet due to capillary and recoil 
forces are investigated for various values of dimensionless evaporation rate, density ratio, 
Reynolds number, disturbance wavenumber, and initial disturbance amplitude. 

Detailed description of the shape evolution of liquid jets with Re=200 and 10 is presented 
in Figs. 2 and 3 for k=0.7,0A and P=0, 10'^, and 3x10'^. The following characteristics for the 
breakup of an evaporating capillary jet can be observed from Figs. 2 and 3: (i) As evaporation 
rate increases the breakup time decreases; (ii) As Reynolds number decreases the breakup time 
increases; (iii) The length and diameter of the liquid ligament decreases with increasing 
evaporation rate and decreasing Reynolds number; (iv) The breakup point is closer to the initial 
trough as evaporation rate increases, due to the decrease in time that nonlinear effects can 
influence jet breakup. A more detailed description of these characteristics is given in the 
following sections. 
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Figure 2; Time evolution of the instability, Re=2QQ, £o=0.5, 5=10’^: (a) k=0.7, P=0; (b) k=0.7, 
p=10'^; (c) k=0.7, p=3xl0'^; (d) k=0.4, p=0; (e) k=0A, P=10'^; (/) k=0.4, p=3xl0'l 



3.1 Location of the breakup point 



Linear theory predicts that the breakup point is always at the trough of the initial 
disturbance where the initial jet radius is minimum. However, Figs. 2 and 3 clearly show that the 
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t=0 




10.05 



11 23 



11.55 



11 65 




1=0 



4.80 



5 09 



5.14 



5 15 



(d) (e) if) 




I 



Figure 3: Time evolution of the instability, /?e=10, 6o=0.5, 5=10’^: (a) k=0.1, P=0; (b) k=0J, 
P=10^ (c) k=0J, p=3xl0'^ (d) k=0A, p=0; (e) k=0A, p=10’^; (/) k=0A, p=3xl0’l 



breakup point moves in the direction of the swell towards the end of the simulation. The 
capillary forces on the jet tend to “squeeze” the liquid in the neck region into the swell. This 
results in the formation of two distinct regions in the jet. The middle region where the curvature 
in the r-z plane is very small, resembling a liquid ligament, and the swell region where the 
curvature in the r-z plane is large. Note from Figs. 2 and 3 that the disturbances stay sinusoidal 
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for a very long time (relative to the breakup time) and agree with linear theory until non-linear 
effects begin to move the location of the minimum radius. 

In order to gain a better understanding of the development of the breakup point of the jet, 
we have plotted the motion of the minimum radius along the jet from time t=0 to the breakup 
time. This plot is shown in Fig. 4 for a jet with /?e=200, £o=0.05, and for three different 

evaporation rates. Two distinct time periods can be identified for each of the curves in these 
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Figure 4: Motion of the minimum radius along the jet axis: (a) P=10'^; (b) p=2xl0'^; (c) P=3xl0'^. 
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figures. The most important being the time that the minimum point stays at z=0. This time 
represents over 90% of the total breakup time in all cases and it first reduces and then increases 
with an increase in the wavenumber. The case with the maximum disturbance growth rate is the 
first to begin movement away from z=0. The second distinct time period represents the rapid 
movement of the minimum jet radius from the neck towards the swell region. This sudden 
relocation of minimum point in the axial direction is larger for smaller wavenumbers. The end 
points on all the curves in Fig. 4 represent the breakup times. There are two interesting features 
to note from Fig. 4. The first is that as evaporation rate increases the location of the minimum 
point is closer to the neck for each wavenumber. This represents the production of a smaller 
middle ligament or satellite drop. The second feature reveals that as evaporation rate increases 
the range of time that each wavenumber breaks away from z=0 reduces. 



3.2 Growth rate of the disturbances 

Linear theory not only defines the region of unstable disturbance wavenumbers but also 
provides their growth rates. These growth rates are useful in estimating the breakup time and 
length. According to linear theory the variation of the logarithmic value of the amplitude of the 
surface disturbances with time is linear. However, for an actual liquid jet this amplitude 
variation is nonlinear at the beginning of the simulation due to the fact that radial velocity of the 
jet starts from zero. The amplitude variation is nonlinear also near the breakup point due to the 
movement of the minimum radius. For a selected few of our simulations we have plotted the 
logarithmic values of the normalized amplitudes of the neck ((Re-rn)/£o), swell ((r^-ReV^o), and 
their difference ((/"s-^nVCo) in Fig. 5. Here, and are the radii of the neck and swell. 
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Figure 5: Variation of the amplitude of the neck (black), swell (red), and their logarithmic 
difference (green), /?e=200, Eo=0.5, s=10'^: (o) k=0.4, p=10 ^ (b) k=0J, p=10 ^ (c) k=0A, 
P=2xl0‘^; (d) k=0.1, (3=2x10 ^ (e) k=0A, (3=3x10'^ (f) k=0.1, p=3xl0 l 



respectively, and Re is the instantaneous radius of the equivalent cylindrical jet, which changes 
during the simulation due to evaporation. Notice in Fig. 5 that the logarithmic variation of 
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amplitude of the neck or swell is not perfectly linear. Figure 5 reveals interesting features for the 
temporal evolution of amplitudes of the neck and the swell with respect to each other. Initially 
the growth rate of the neck point is larger than the growth rate of the swell point, since a radial 
displacement in the neck region corresponds to a smaller radial displacement in the swell region 
for the same volume displacement. Due to the formation of a liquid ligament, i.e. minimum 
radius point moving away from the neck, the growth rate of the amplitude of the neck decreases 
eventually dropping below that of the swell for the longer wavelengths. This is clearly shown in 
Fig. 5a,c,e, which also indicates that as we increase the evaporation rate the time after the cross 
point is reduced. This represents the formation of a smaller liquid ligament, which can be seen 
in Figs. 2 and 3. Figure 5b,d,f is representative of smaller disturbance wavelengths. Due to the 
reduced time of breakup non-linear effects are only visible in the last instances of instability. 
Comparing Fig. 5a,c,e with 5b,d,f we see that the shape of the difference curve varies only in 
slope when a change in wavenumber is made. 

Similar to previous research in this area we use the linear region of the difference curve 
to provide values that can represent the growth rate for each wavenumber. Therefore, a line is fit 
to the difference between the amplitudes of the neck and swell and its slope is measured. The 
line is fit only to the middle region of each curve and the end portions are ignored. The 
dispersion curves are obtained using the calculated values of the growth rate for different 
wavenumbers and for various values of p. Re, 8o, and s. The data are plotted in Fig. 6. As 
mentioned earlier Fig. 6a shows quantitatively that as we increase the evaporation rate the 
growth rate increases therefore reducing the time it takes for the jet to breakup into droplets. All 
curves in Fig. 6a give the maximum growth rate at k=0J, similar to linear theory. Figure 6b 
presents the effects of Reynolds number on growth rate of disturbances. Here we can clearly 
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Figure 6; Growth rate of disturbances: (a) /?e=200,eo=0.05,5=10'^; {b) £o=0.05,P=3xl0'^,5=10‘^; 
(c) ^c=200,P=3x10'^5=10'^; {d) /?e=200,eo=0.05,p=3xl0'^; {e) Lian & Reitz (1993). 



identify for evaporating jets that as we decrease the Reynolds number the growth rate of 
disturbances decreases due to the more viscous nature of the jet not allowing the fluid to be 
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displaced as easily. Also from this graph we see a movement of the maximum growth to a lower 
wavenumber for the case of highly viscous fluids. Figure 6c demonstrates the effects of initial 
disturbance amplitude on growth rates. As mentioned earlier, linear theory assumes infinitesimal 
initial disturbance amplitude. Therefore, for large initial disturbance amplitudes the theoretical 
solutions result in large errors. We see from Fig. 6c that as we increase this value growth rate 
decreases at each wavenumber. However, this does not mean that the jet will take longer to 
breakup for large initial disturbance amplitudes. This figure only gives meaning to the rate and 
not the overall distance that the radius has to decrease in order to pinch off the jet. Figure 6d 
provides insight into the effect of the density ratio on growth rate. Here we clearly see that as we 
increase the density of the gas or decrease the density of the liquid (i.e. larger s) the growth rate 
of the evaporating jet decreases. This can be explained through the reduction of the recoil force, 
where the net momentum imparted on the surface is diminished by increasing s (see (13)). 

Figure 6e shows the results from Lian and Reitz’s (1993) linear stability analysis for an 
evaporating viscous jet with 5=0.25x10’*, Weg=10’^, and Z=10’'*, where We^ represents the 
Weber number and Z represents the Ohnesorge number. Since we do not solve for the gas phase 
we must convert their formulation to agree with our input parameters. For the case shown in Fig. 
6e this would result in a Re=20,000 due mainly to their very small value for density ratio (e.g. s 
for air and water would be approximately 1.2x10’ ). Using this Reynolds number and similar 
values of s and (3 we attempted to correlate a direct comparison but was unable to reproduce their 
results. The calculated growth rates in our simulations were very similar to cases with no 
evaporation. Only when we increased P to around 5x10’^ did we obtain growth rates larger than 
that of non-evaporating cases. Although we could not reproduce the exact results from the 
analytical study, even using the same parameters, we both conclude the same general trends for 
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growth rate. As evaporation rate increases growth rate increases for every wavenumber. One 
possible reason for the discrepancies between our results and theirs is the linear theory’s 
assumption of infinitesimal disturbance amplitude. For non-evaporating jets the growth rate 
does not vary up to an initial disturbance amplitude of approximately 0.05. However, Fig. 6c 
clearly shows that the growth rate varies in this region for evaporating jets. From this one may 
question the validity of the normal mode method for evaporating jets. Another possible reason 
for the differences can be tracked to their initial conditions where at t=0 the free surface is 
moving but there is no internal flow field. We could not produce such initial condition in our 
numerical simulations. 



3.3 Breakup time and volume 

Linear theory predicts that the minimum breakup time occurs at k=0.691, which 
corresponds to the wavenumber with the maximum growth rate. However, this is only true for 
£o— >0, otherwise the minimum breakup time occurs at higher values of k as explained by 
Chaudhary and Redekopp (1980). In the simulations presented in Fig. 7a,b,c,d the initial 
disturbance amplitude equals 0.05 therefore most of our curves give a minimum breakup time for 
k=0.8. From Fig. 7a we can quantitatively deduce that as we increase the evaporation rate we 
decrease the breakup time, which has been demonstrated by other graphs presented earlier. This 
is due to the fact that the addition of recoil force has a destabilizing effect for low speed jets, 
which we are considering in this study. Also evaporation decreases the mean radius of our jet 
with time. Figure 7b clearly shows that as we decrease the Reynolds number breakup time 
increases, which is due to the additional viscous damping. This result is also in agreement with 
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Figure 7: Breakup time of jet: (a) /?e=200,Eo=0.05,5=10'^; (b) eo=0.05,P=3xl0'^,5=10'^; (c) 
/?e=200,(3=3xl0 ^5=10 ^ (d) /?e=200,Eo=0.05,p=3xl0‘l 



previous observation for non-evaporating jets in the Reyleigh regime (Mashayek and Ashgriz 
1995). Another important point to note is that the minimum breakup time shifts to lower 
wavenumbers as Re is significantly reduced. This is in accordance to linear theory which states 
that the maximum growth rate shifts toward the lower wavenumbers for more viscous liquid jets. 
Even though our cases involve evaporation. Fig. 7c still re-emphasizes Chaudhary’s results by 
showing that as Eo increases from zero the minimum breakup time occurs at higher 
wavenumbers. For example, for Eo=0.5 the minimum breakup time occurs at k=0.9 but for 
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£o=0.001 this happens at ^=0.8. Clearly, from Fig. 7c we notice as initial disturbance amplitude 
increases breakup time decreases. From inspection of Fig. 7d we can dissimulate that an 
increase in density ratio increases breakup time due to the diminishing of the recoil force. 



Breakup volume, the instantaneous combined volume of both the main drop and the 
satellite drop at breakup, is primarily dependent on jet breakup time. If it takes longer for jet 
breakup to occur then the jet is allowed more time to evaporate. Beginning with an initial 
normalized jet volume of one. Fig. 8 shows the final volume of the jet while varying different 
parameters. Analysis of the results from each graph gives similar conclusions to that described 
for the breakup time. 
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Figure 8: Breakup volume of jet: {a) /?e=200,£o=0.05,5=10'^; ib) £o=0.05,P=3xl0'^,5=l0'^; (c) 
/?^=200,p=3xl0 ^5=10'^ {d) /?e=200,£o=0.05,p=3xl0'l 



22 






}*' .r U'l.Xij 






tir 



•if n; 



■ -i. , 



t 



r 



! 



s 






'•t M»u* * ri»faw?lt 






i ii.. 



A fi: Jfr- 






3.4 Drop size at breakup 



The control of satellite drops has been one of the primary driving forces behind nonlinear 
jet instability studies. For non-evaporating jets the total volume of a single wavelength of fluid 
is constant throughout the breakup process and equal to the combination of the main drop 
volume and satellite drop volume at breakup. However for evaporating jets the volume changes 
in time and the total volume of the jet at any instant is equal to the initial jet volume minus the 
total amount of evaporation up to that point. At breakup, the volume of an evaporating jet is 
equal to the combined volume of both the main and satellite drops, where this total could be 
significantly less than the initial volume depending on the evaporation rate and the breakup time. 




Figure 9: Variation of drop size at breakup, /?e= 200 ,£o= 0 . 05 , 5 = 10 '^. 
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Figure 9 plots the main and satellite drop radii versus wavenumber for p=10' ,2x10' ,3x10' . 
From this figure we see that, for small evaporation rate, as wavelength increases the drop size 
generally increases, which is similar to non-evaporating jets. Flowever, large evaporation rates 
and small wavenumbers (i.e. k=0.3 for (3=2x10'^ and ^=0.3, 0.4 for (3=3x10'^) we notice a 
decrease in both main and satellite drop size due to the increased breakup time, which allows 
more vapor to leave the surface and reduce the overall jet volume. Two general conclusions 
from Fig. 9 are; (i) as evaporation rate increases drop size decreases for all disturbance 
wavenumbers, and (ii) comparing the percent reduction in drop size evaporation affects the 
satellite drop more than the main drop (e.g. by increasing (3 by a factor of 3 at k=0.3 the main 
drop radius reduces 38% and the satellite drop radius reduces 50%, or at A:=0.9 this reduces the 
main drop by 10% and the satellite drop by 29%). 



3.5 Evaporation flux on free surface 

In our evaporation model we approximate the unit area evaporation rate of our perturbed 
jet surface with that of a sphere with the same local radius of curvature. Curvature of the jet 
surface is determined by the height function using (15). To be in agreement with (9) the 
curvature of the surface was monitored throughout simulations to prohibit negative values in 
evaporation cases. Figure 10 shows the curvature of the jet versus the axial coordinate z for the 
simulations presented in Fig. 2. In every case the curvature increases in the neck region as we 
approach the breakup time. This translates to an increase in the evaporation mass flux as shown 
in (26). That is, the evaporation flux at the neck is higher than that at the swell, leading to 
further pinching of the neck. Also the recoil force is stronger in the neck region than that at the 
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Figure 10: Variation of the curvature of the jet, ^e=200, £o=0.5, 5=10'^: (a) k=0.7, P=0; (b) 
k=0.4, p=0; (c) k=0.7, P=lxl0'^ (d) k=0.4, p=lxl0'^; (e) k=0.7, p=3xl0'^; (/) k=0A, p=3xl0 l 



crest, causing the liquid to be squeezed into the crests and thus accelerating the growth of the 
disturbance at the neck. The point of highest curvature in each set of data represents the 
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instantaneous minimum radius of curvature. At the breakup time we see a drastic increase in the 
curvature which is associated with the breakup point. By investigating Fig. 10a,c,e we can easily 
see that the breakup point is closer to the neck for increasing evaporation rates. 



3.6 Pressure and velocity distributions 

To provide further insight into the dynamics of the jet the velocity and pressure 
distributions within the liquid jet have been investigated. Figure 11 shows the pressure 
distribution and velocity vectors at three different non-dimensional times for a non-evaporating 
jet. At all the times shown we can easily see the flow of the fluid from the neck to the swell. 
With increasing surface deformation a significant increase of velocity in the neck region takes 
place. Initially the pressure distribution is uniform within the jet but as the velocity field 
develops we can clearly identify areas of low and high pressure. At the final time shown, which 
is near the breakup time, we can distinguish a high-pressure region at the top of the swell due to 
the conversion of velocity to pressure. The middle ligament exhibits a lower pressure where the 
velocities are higher. For comparison Fig. 12 gives the results for an evaporating jet with similar 
parameters as Fig. 11. From the legend we see that the evaporating case produces a higher 
overall pressure, which is due to the contribution of both capillary and recoil forces. We can 
clearly identify the recoil force in both the middle ligament and the main drop of our jet where 
we see a higher pressure along the surface. It should be mentioned that in Figs. 1 1 and 12 the 
velocity vectors have been scaled in each of the 3 non-dimensional times for clarity. 

By increasing the disturbance wavelength and decreasing the initial disturbance 
amplitude our code has been able to predict multiple satellite drops, which has been investigated 
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Figure 1 1 : Pressure distribution and velocity vectors for a non-evaporating jet, Re=200, k=0J, 
£o=0.05, P=0, 5=10 l 



experimentally by Vassallo and Ashgriz (1991). Figures 13 and 14 give insight into the 
differences seen between non-evaporating and evaporating jets under these conditions. From the 
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Figure 12: Pressure distribution and velocity vectors for an evaporating jet, Re=200, k=0.7, 
£o= 0.05, P=3xl0■^ ^=10 l 



non-evaporating case. Fig. 13, we notice that nonlinear effects have drastically influenced the 
surface shape near breakup. Here we can identify a main drop and, possibly, multiple satellite 
drops which after breakup would either combine and oscillate or separate. Investigating Fig. 14 
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Figure 13: Formation of multiple satellite drops for a non-evaporating jet, /?e=200, k=0.3, 
eo=0.001,p=0,5=10^ 




Figure 14: Formation of multiple satellite drops for an evaporating jet, Re=200, k=03, £o=0.001, 
p=10 ^ 5=10 ^ 



for evaporating case, we see a different surface shape due to the relocation of the breakup point 
close to the initial neck. This jet creates one small liquid ligament and two larger bulged 
volumes of liquid, which would continue to interact after breakup. 



4. Conclusion 

Numerical simulation is used to investigate the instability of evaporating capillary jets. 
The evaporation model is based on the rf^-law, widely used for spherical drops. To implement 
this model for perturbed jets, the conventional (f-law is modified by substituting the local radius 
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of the surface for the radius of a sphere. The need for a solution for the ambient gas is 
eliminated considering small density ratios. Numerical simulations are carried out using a 
Galerkin finite element method with penalty function. The motion of the free surface is tracked 
using the Height-Flux Method developed by Mashayek and Ashgriz (1993). With this technique, 
the evaporation at the surface is easily handled by modifying the amount of fluid within each 
subvolume. Preliminary simulations were performed to establish the required mesh size and 
time increment for various cases. 

In non-dimensional form, the formulation involves five parameters; the amplitude of the 
initial surface disturbance, the wavenumber of the initial surface disturbance, the Reynolds 
number, the non-dimensional evaporation rate, and the density ratio. A wide range of variation 
for each of these parameters is considered in numerical simulation. The results were discussed 
for various cases leading to the foremost conclusion that evaporation clearly has an effect on jet 
instability, mainly increasing the growth rate of the disturbances by increasing the evaporation 
flux at the neck. Thus, evaporation causes the jet to break up sooner and reduces the size of both 
main and satellite drops. An increase in internal pressure is also noticed with evaporating cases 
due to the contribution of both capillary and recoil forces. Due to the limited number of 
evaporating jet instability studies the majority of our conclusions could not be compared to 
experimental or analytical results. 
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